library(hydroGOF)
library(clusterSim)

data=read.csv("C:/data/mod_merge3.csv")
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")
for (i in 2010:2014) {
	testyear = i
	trainData=subset(data,data$year!=testyear)
	testData=subset(data,data$year==testyear)
	data1=trainData
	
	model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
	pm25model=predict(model1,testData)
	
	
	model_re=sum(abs(testData$pm25-pm25model)/testData$pm25)/nrow(testData)*100
	model_rmse=rmse(testData$pm25,pm25model)
	model_r=cor(testData$pm25,pm25model)
	
	print(paste("Year model 1: ",testyear))
	print(paste("Number of training: ",nrow(data1)))
	print(paste("Number of testing: ",nrow(testData)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))
	print("------------------")
	
	nSample=nrow(data1)
	nPredict= length(model1$coefficients)-1
	
	cutoff=4/(nSample-nPredict-1)
	cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
	cookValue=cooks.distance(model1)
	data1["cookValue"] = cookValue
	
	data2=subset(data1,data1$cookValue<=cutoff)
	data3=subset(data1,data1$cookValue<=cook_thresold)
	
	model2=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data2)
	pm25model=predict(model2,testData)
	
	
	model_re=sum(abs(testData$pm25-pm25model)/testData$pm25)/nrow(testData)*100
	model_rmse=rmse(testData$pm25,pm25model)
	model_r=cor(testData$pm25,pm25model)
	
	print(paste("Year model 2: ",testyear))
	print(paste("Number of training: ",nrow(data2)))
	print(paste("Number of testing: ",nrow(testData)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))
	print("------------------")
	
	model3=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data3)
	pm25model=predict(model3,testData)
	
	
	model_re=sum(abs(testData$pm25-pm25model)/testData$pm25)/nrow(testData)*100
	model_rmse=rmse(testData$pm25,pm25model)
	model_r=cor(testData$pm25,pm25model)
	
	print(paste("Year model 3: ",testyear))
	print(paste("Number of training: ",nrow(data3)))
	print(paste("Number of testing: ",nrow(testData)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))
	print("------------------")
}

